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IMPLEMENTATION  OF  A  DISCONTINUOUS  GALERKIN  DISCRETIZATION  OF 
THE  CONSERVATION  OF  MASS  EQUATION  IN  QUODDY 

1.  INTRODUCTION 

In  this  effort,  a  discontinuous  Galerkin  (DG)  formulation  is  developed  for  the  primitive  2-D  unsteady 
depth-averaged  conservation  of  mass  equation.  The  formulation  is  implemented  into  QUODDY5, 
developed  by  Lynch  and  Werner  [1]  at  Dartmouth  University  and  capable  of  solving  the  3-D  shallow 
water  hydrodynamics  and  temperature  and  salinity  transport. 

Traditional  continuous  Galerkin  finite  elements  have  historically  not  performed  well  when  applied  to 
the  primitive  form  of  the  shallow  water  conservation  of  mass  equation  due  to  spurious  oscillations 
introduced  into  the  solution.  In  an  effort  to  eliminate  the  spurious  oscillations.  Lynch  and  Gray  [2] 
developed  the  generalized  wave  continuity  equation  (GWCE).  The  method  eliminated  the  spurious 
oscillations  and  has  been  adopted  in  several  shalow  water  hydrodynamic  codes,  for  example  QUODDY 
[1]  and  ADCIRC  [3],  The  term  discontinuous  Galerkin  (DG)  method  refers  to  several  finite  element 
methods  that  use  discontinuous  discrete  approximating  spaces,  such  as  the  Bassi  and  Rebay  method  [4], 
the  Local  Discontinuous  Galerkin  (LDG)  [5]  method,  the  Oden,  Babuska,  and  Baumann  method  [6],  the 
interior  penalty  Galerkin  methods  of  Wheeler  [7],  Douglas  and  Dupont  [8],  and  the  NIPG  methods  [9,10]. 
Cockbum  et  al.  performed  an  extensive  study  of  the  DG  method  for  hyperbolic  conservation  laws  in  a 
series  of  papers  [11-15].  Arnold  et  al.  [16]  present  a  general  framework  of  these  methods.  Application  of 
these  methods  to  a  wide  variety  of  problems  can  be  found  in  Ref.  17. 

A  unique  property  of  the  DG  method  is  that  the  DG  solution  satisfies  the  mass  conservation  equation 
locally  element  by  element.  In  addition,  the  flexibility  of  DG  methods  facilitates  dynamic  h-p  adaptivity. 
The  basis  functions  are  defined  as  piecewise  discontinuous  polynomials  on  each  element,  so  the  order  of 
polynomials  can  be  refined  locally  (p-adaptivity).  Additionally,  elements  can  be  nonconforming,  which 
facilitates  local  mesh  refinement  (h-adaptivity).  These  features  can  be  exploited  separately  or  in 
combination  (h-p  adaptivity)  to  locally  refine  the  approximation  in  regions  of  high  gradients.  Moreover, 
because  they  are  highly  local,  DG  methods  only  require  communication  between  elements  that  share 
faces;  thus  they  are  well-suited  for  parallel  computation. 

2.  ORIGINAL  FINITE  ELEMENT  FORMULATION  IN  QUODDY 

QUODDY  is  capable  of  computing  3-D  shallow  water  flow,  as  well  as  transport  of  temperature  and 
salinity.  The  system  of  equations  for  computing  the  velocity  field  implemented  in  QUODDY  are  the 
GWCE,  the  3-D  primitive  momentum  equation,  and  3-D  incompressible  conservation  of  mass  equation. 
The  GWCE  is  used  to  compute  the  elevation  given  the  depth-averaged  velocity  components.  The  3-D 
primitive  momentum  equations  are  used  to  compute  the  3-D  horizontal  components  of  velocity,  and  the 
3-D  incompressible  conservation  of  mass  equation  is  used  to  compute  the  vertical  velocity.  Reference  1 
provides  a  detailed  discussion  of  the  governing  equations  and  finite  element  method  implementation. 
QUODDY  is  capable  of  incorporating  land,  velocity,  and  elevation  boundary  conditions  into  the  solution. 
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The  numerical  procedure  currently  in  QUODDY  is  as  follows: 

1 .  Given  an  initial  depth  H  and  depth-averaged  velocity  field  (U,  V),  update  the  elevation  by  solving 
the  GWCE  using  a  three-level  time  differencing.  The  solution  is  semi-implicit  in  that  all  of  the 
linear  terms  are  solved  implicitly,  but  all  of  the  nonlinear  terms  involving  velocity  are  lagged. 

2.  Using  the  new  computed  depth,  compute  the  components  of  velocity  for  the  entire  3-D  domain 
using  a  two-level  time  differencing. 

3.  Compute  the  component  of  velocity  such  that  the  velocity  field  is  divergence  free,  i.e.,  V  *  V  =  0. 

3.  DISCONTINUOUS  GALERKIN  FORMULATION 


In  this  section,  a  DG  formulation  on  triangular  elements  is  developed  and  presented.  The  method 
begins  by  forming  a  weighted  residual  statement  and  integrating  over  a  typical  element.  Summing  over  all 
elements  produces  a  system  of  equations  to  advance  the  solution  in  time.  Integrating  the  test  and  basis 
functions  over  each  element  produces  volume  and  interface  integrals  that  must  be  evaluated.  In  this 
formulation,  the  volume  integrals  are  evaluated  analytically  and  the  element  interface  integrals  are 
evaluated  using  Gaussian  quadrature. 

3.1  Problem  Statement 


The  primitive  form  of  the  2-D  depth- averaged  inviscid  shallow 
conservative  form  as 


water  equations  is  written  in 


where 
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are  the  3x1  vectors  that  represent  conservation  of  mass,  x- momentum,  and  v-momentum,  respectively. 

This  effort  focuses  on  incoiporating  the  conservation  of  mass  into  QUODDY  using  the  DG  method. 
However,  the  continuous  Galerkin  method  already  implemented  into  QUODDY  is  retained  for 
conservation  of  momentum. 

Before  developing  the  numerical  formulation,  element  nomenclature  used  throughout  the  remainder  of 
this  work  is  defined.  Each  element  contains  three  nodes  and  three  edges  as  shown  in  Fig.  1(a).  The  edges 
between  each  element  are  denoted  by  Y k .  Each  interior  edge  is  attached  to  two  elements,  defined  as  the 

left  ( L )  and  right  (R)  element.  The  unit  normal  to  an  edge  nk  is  defined  such  that  it  always  points  outward 

from  the  left  element.  This  is  described  in  Fig.  1(b).  Edges  on  the  boundary  are  only  attached  to  the  left 
element.  The  boundary  condition  replaces  the  right  element. 
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Fig.  1  — (a)  Element  and  (b)  edge  nomenclature  for  typical  interior  elements 


The  DG  method  involves  jumps  and  averages  across  edges.  The  jump  and  average  of  any  quantity  (a) 
across  edge  k  are  defined,  respectively,  as 


[a\k  =  (aR  -  aL ) 
{a)k  =- {aL  +cir)- 


Multiplying  conservation  of  mass  by  a  test  function  w  and  integrating  over  a  typical  element 
produces  the  semidiscrete  weak  statement  of  conservation  of  mass  and  is  written  as 


d_ 

dt 


j  W:HjdQ.j  + 


Q, 


3 

I 


k=\ 


J  w(F-nkfiyTk-  j  (F-Vw)dnj=  0, 
r  k  cij 


(4) 


where  T  k  are  the  edges  of  element  j,  and  F  =  (UI I )  i  +  ( VII  j  j  is  the  flux  vector.  In  Eq.  (4),  the  normal 

nk  0  =  nj  +  nvj  is  the  unit  outward  normal  vector  to  element  j.  Therefore,  F  -nko  is  the  outward  flux 

normal  to  the  element  edge  k.  As  noted  previously,  normals  in  the  implementation  of  the  algorithm, 
defined  as  per  Fig.  1(b),  point  from  L  to  R.  Therefore,  for  all  right  elements,  the  normal  for  an  edge  must 
be  multiplied  by  minus  one.  This  prevents  the  need  to  define  two  normals  for  each  edge. 

3.2  Basis  and  Test  Functions 

In  the  current  work,  linear  basis  and  test  functions  are  used.  Linear  basis  and  test  functions  are  defined 
for  the  /h  element.  Linear  basis  functions  provide  three  degrees  of  freedom  per  element.  Several 
formulations  are  possible  depending  on  the  degrees  of  freedom  chosen.  The  general  forms  of  the  basis 
and  test  functions  are 

Hj  (*)  =  Z  <H <t>*  wj(x)=  Z  bi (5) 

i= 1  i=l 

The  specific  form  of  the  basis  functions  (j)/  depends  on  the  degrees  of  freedom  chosen.  Commonly,  the 
degrees  of  freedom  are  the  element  average  and  gradient,  the  three  element  vertices,  or  the  three  element 
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edge  midpoints.  In  this  work,  two  formulations  are  investigated:  element  average  and  gradient  (denoted  as 
Method  1)  and  element  vertices  (denoted  as  Method  2). 

The  functions  4*/  and  coefficients  a,  are  given  in  Table  1.  The  coefficients  b,  (Eq.  (5))  are  arbitrary 
numbers  that  are  eliminated  from  the  final  set  of  equations.  The  element  centroid  is  written  as  ( x ,  y ) . 


Table  1  — Test  Functions  and  Coefficients 


i 
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and  Gradient 
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a,- 
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at 
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3  2Ae 
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dir 
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dx 

3  2Ae 
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h3 

1  !  {y\-yi)(x-x)(xx-X2)(y-y) 

3  2Ae 

3.3  Linear  Integration  Rules  for  Triangular  Elements 

Applying  the  linear  basis  and  test  functions  of  the  DG  method  to  triangular  elements  produces 
integrals  of  the  form 


Ix=\{x-x)dA  Iy  =  J(  v- v)  dA  Ixy  =  \(x-x){y-y)dA. 

AAA 

The  integrals  are  evaluated  as 


IX  =  A 


ly=A 


Jxy=A 


n^+x2+xl)-^x2  ^ 

T(ti2  +y2+yl)~\y2 


— (*lTl  +  *2k2  +  x3T3  )  -  ~xy 


(6) 


(7) 


3.4  Edge  Fluxes 

The  weak  statement  of  conservation  of  mass  includes  integrals  of  the  fluxes  at  the  edges  of  each 
element.  Upwinding  is  achieved  by  replacing  the  normal  flux  in  the  edge  integral  of  Eq.  (4)  with  an 
upwinded  numerical  flux  function  normal  to  the  element  edge,  denoted  by  fn.  The  numerical  flux  function 
is  computed  by  solving  a  Riemann  problem  normal  to  each  element  interface.  With  this  substitution,  Eq. 
(4)  becomes 
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{  WjHjdD.j =  }  (F-Vw)dD.j-Z  1  wfndYk. 

fi7.  n J  k=i  r, 


This  can  be  written  in  matrix  form  as 


[M] 


d[  H] 
dt 


[V][H]-[F], 


(8) 


(9) 


where  [M]  is  the  mass  matrix,  [V]  is  a  3  x  1  matrix  involving  the  volumetric  flux  integral,  and  [F]  is  a 
3  x  1  vector  involving  the  edge  flux  integral.  [H]  is  a  3  x  1  vector  involving  the  degrees  of  freedom  for 
each  element.  The  specific  form  of  the  vectors  and  matrices  depends  on  the  method  chosen. 

3. 4. 1  Riemann  Problem 


The  jumps  in  the  dependent  variables  across  element  interfaces  produce  a  Riemann  problem  normal  to 
each  element  interface.  The  governing  equations,  as  represented  by  Eqs.  (1)  and  (2),  are  linearized  by 
Roe-averaging,  and  the  solution  to  the  linearized  Riemann  problem  computed.  Once  the  solution  to  each 
local  Riemann  problem  is  obtained,  the  edge  fluxes  are  computed.  The  solution  to  the  linearized  Riemann 
problem  is  well  known  (see,  for  example,  Leveque  [18]). 

The  solution  to  the  normal  Riemann  problem  is  written  in  terms  of  the  eigenvalues  and  eigenvectors  of 
the  normal  Jacobian  matrix  of  Eq.  (1),  which  is  written  in  quasilinear  form  as 


dQ 

dt 


+  [A] 


7T+[5] 

ox 


^-  =  S(x,t), 
oy 


(10) 


where  \A\  and  [B]  are  3x3  matrices. 


The  Roe-averaged  normal  Jacobian,  [A]„  =  [A]nx  +  [B]ny,  is  written  as 


0 


[Aln  = 


(-U2nx  +  gHnx-UVny) 
(-V2ny  +  gHny-UVnx) 


(2  Unx  +  Vny)  ( Uny ) 
(Vnx)  (Unx  +  2Vny) 


The  eigenvalues  of  \A]„  are 


^1  =  un  ~  sfgH 

X2  =  un 

h  =  K  + 


(11) 


(12) 


where  un  -  Unx  +  Vn  is  the  velocity  normal  to  the  edge  and  the  overbars  denote  Roe-averaged  values. 
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The  corresponding  right  eigenvectors  are 
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U  +  s]gHnx 
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ny 
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The  Roe-averaged  values  are  defined  for  an  edge  k  as 

_  Hl  +  Hr 
2 

-  unL  ^h  +  unR 


(13) 


(14) 


The  subscript  (L,R)  denotes  computing  the  given  quantity  within  the  left  or  right  element  at  edge  k. 
The  value  of  the  numerical  flux  function  for  conservation  of  mass  on  given  edge  depends  on  the  signs  of 
the  eigenvalues  at  the  edge.  The  numerical  flux  function  for  conservation  of  mass  is  written  as 


fn 


fn  ~  unL^ L 
fn  =  unRHR 

X  — 

unLHL  +  7T  _*  7T  [^3 (HR  ~Hl)-(uhRHR  ~uhLHl)~\ 


>  0  and  X3  >  0 
kj  <  0  and  X3  <  0 
<  0  and  X3  >  0 


|  Critical  and  Supercritical 


(15) 


}  Subcritical 


4.  SOLUTION  ALGORITHM 


The  solution  to  Eq.  (8)  involves  evaluating  the  integrals  using  the  analytical  integration  rules  presented 
in  Section  3.3  for  the  volume  integrals,  and  Gaussian  quadrature  for  the  edge  integrals.  The  time 
derivative  is  approximated  using  a  second-order  Runge-Kutta  method  between  time  levels  n  and  n+  1 . 

4.1  Temporal  Discretization 

Two  variations  of  a  second-order  Runge-Kutta  method  were  investigated.  Each  is  a  predictor-corrector 
type  discretization.  The  first  method  advances  the  solution  from  time  level  to  time  level  as  follows: 

Hn+ji  =  Hn  +  ^-L(Hnyn,Fn) 


H"+1 


=  H"  +  AtL\ 


«+-'■ 

H  2 


V  2 


,F  2 


v 


j 


(16) 


The  second  method,  known  as  Henn’s  method,  advances  the  solution  through  the  following  predictor- 
corrector  sequence: 


Implementation  of  a  Discontinuous  Galerkin  Discretization 


7 


H*=H"+AtE(H",V",F") 

H"+l  =  H"  +  y  (L(  H" ,  V"  ,F”)  +  L(  H*  ,V*,F*))’ 

where  f(H,V,F)  are  written  as 

Z,(H,V,F)  =  M_1  VH  -  — M  1  F. 

A 

As  discussed  below,  the  second  method  is  somewhat  more  stable  than  the  first. 

A  3  x  3  system  of  algebraic  equations  is  produced  for  each  element  to  update  the  value  of  Hj.  Note  that 
each  element  is  coupled  only  to  its  neighboring  elements  via  the  edge  flux  integrals. 

4.2  Volume  Integrals 

The  volume  integral  involving  the  time  derivative  at  time  produces,  for  methods  1  and  2,  respectively, 
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Similarly,  the  volume  integral  of  the  flux  produces, 

0 


J  (. Fn-Vw)dQ 


Q, 


=  AjV"H  =  Aj 


0 


U'j 0 


' 4  8UJ  +  4y  duf 
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A  dx  A  dy 


j  0 


I  dV '■  I  dVn 

Jx  uv  J  +  xy  ur  j 


A  dx  A  dy 
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Ixy  dUj  +Iy  dU) 


A  dx  A  dy 
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xy  uv j  Ay  uv j 
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HJ0 

dHj 

\ 

dx 

dHj 

)  _ 

dx 

(y2 -y3)(  2U'l  +U2+U%)  (y2 -y3)(U{1  +2UJ{  +  U%) 
(y3 -yi)(2 U?+U%  +U,{)  (_y3  ~y\)(U\  +  2U2  +  U3) 
(Tl -y2)(2U?+UZ  +C/3b)  Cyj -y2Wx  +2U%  +U%) 


A  dx  A  dy 


(y2  —  y3)(U\  +2u2  +u3) 
(y3  -yi)(U\  +2u2  +u3) 
(ti  ~  y2)  (JJi  +  2U2  +  u3 ) 


(20) 


~H{ 

h2 

H3_ 

In  Eq.  (20),  the  velocities  have  been  projected  from  nodal  degrees  of  freedom  to  DG  degrees  of 
freedom.  The  projection  algorithm  is  described  in  Section  5.2. 
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4.3  Interior  and  Boundary  Edge  Integrals 

For  a  given  edge  k  attached  to  element  j,  the  edge  flux  integral  is  written  as 

3  K 

F-S  <l>2  nTj.  (21) 

k= i  r,  i 

4  uhj 

The  edge  flux  integrals  are  evaluated  using  Gaussian  quadrature.  Both  methods  produce  integrals  of 
the  following  form: 

1  ni  1 

(x-xj)  dTk  =  }  /"  x  d Tk-  xj  }  f„dTk.  (22) 

(y-yj)  \  Tk  [y\  [yj} 

Note  that,  in  general,  the  integrand  depends  not  only  on  the  value  of  the  edge  flux,  but  also  on  the 
values  [  x  -  T/  j  and  (  y  -  T/  j  for  the  given  element  j.  Thus,  the  same  edge  integral  will  have  two  different 

values  for  the  two  elements  sharing  the  edge.  For  computational  efficiency,  the  integrals  given  by  the  first 
expression  on  the  RHS  of  Eq.  (22)  are  evaluated  using  Gaussian  quadrature  for  each  edge,  and  the  total 
integral  for  an  edge  connected  to  element  j  is  found  by  summing  the  two  terms  on  the  RFIS  of  Eq.  (22). 
Using  Gaussian  quadrature,  this  is  evaluated  as 


i  N  n,rn 

I  fn  X  dTk=^-  twm  xmf”n 
4  LtJ  ymf-n 


where  Np  is  the  total  number  of  integration  points,  (xm,  ym)  are  the  Gaussian  integration  points,  wm 
represents  the  Gaussian  weighting  factor,  and  Lk  is  the  length  of  the  edge.  This  implementation  uses  a 
two-point  Gaussian  quadrature. 

For  the  interior  edges,  the  flux  is  evaluated  using  the  Roe-averaged  variables  based  on  the  left  and 
right  states  adjacent  to  the  edge  and  substituting  these  into  Eq.  (15).  The  boundary  edges  are  evaluated 
similarly,  but  the  right  state  is  replaced  by  an  appropriate  boundary  condition. 

There  are  five  possible  boundary  conditions  implemented  into  QUODDY  implementation  of  the 
boundary  condition  for  each  case  described  below. 

1.  Land:  Zero  normal  flow  boundary  condition  is  applied  to  external  boundaries. 

2.  Island:  Zero  normal  flow  boundary  condition  is  applied  to  internal  boundaries. 

Land  and  island  boundaries  are  zero  normal  flow  boundaries.  At  these  boundaries,  the  Roe-averaged 
normal  velocity  is  set  to  zero.  This  is  accomplished  by  setting  the  right  value  of  the  normal  velocity  equal 
to  the  negative  of  the  left  value  and  setting  the  right  value  of  the  elevation  equal  to  the  left  value,  i.e., 
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UnR  llnl 
Hr=Hl 


(24) 


at  land  and  island  boundaries. 

3.  Non-zero  normal  velocity:  Value  of  velocity  is  specified  at  given  boundary. 

At  non-zero  normal  velocities,  the  boundary  specified  normal  velocity  replaces  the  right  value  as 

unR  =  un,BC  •  (25) 

4.  Geostrophic  outflow:  Elevation  is  specified  at  the  given  boundary. 

5.  Surface  elevation:  Elevation  is  specified  at  the  given  boundary. 

Since  geostrophic  outflow  and  surface  elevation  boundary  conditions  both  specify  elevation,  they  are 
handled  the  same  way  in  the  DG  formulation.  The  specified  surface  elevation  at  the  boundary  replaces  the 
right  value  on  elevation  prescribed  edges  as 


hr  ~hbc- 

The  velocity  on  the  boundaiy  is  found  by  a  conservation  of  mass  constraint,  such  that 


unR  ~ 


'jh' 

\HBC  J 


ln,BC  ■ 


(26) 


(27) 


4,4  CFL  Condition 

The  explicit  formulation  of  the  DG  method  restricts  the  time  step  via  a  Courant-Friedrichs-Levy  (CFL) 
condition.  The  courant  number  is  defined  for  each  edge  k  and  the  time  step  is  limited  by  the  maximum 
courant  number  in  the  mesh.  The  courant  number  is  defined  as 


At  |  A,  |  k  ,max  <  j 

min  (Al,Ar) 


VTZ-  gQ, 


where 


A  is  the  maximum  absolute  value  of  At  i  =  1, 2, 3  on  each  edge. 

A:,  max 


(28) 


4.5  Slope  Limiting 

The  stability  of  the  DG  method  cannot  be  ensured  without  incoiporating  some  type  of  slope  limiting 
into  the  algorithm.  Two-dimensional  slope  limiters  are  a  subject  of  current  research.  Cockbum  and  Shu 
[14]  develop  a  total  variation-bounded  (TVB)  slope  limiter  valid  for  linear  approximations  on  rectangles 
and  triangles  with  certain  restrictions.  Floteit  et  al.  [19]  propose  a  new  slope  limiter  they  claim  achieves 
limiting  without  introducing  excessive  numerical  viscosity.  In  that  work,  they  also  review  the  slope 
limiter  of  Cockbum  and  Shu.  In  this  effort,  the  slope  limiter  of  Cockbum  and  Shu  is  currently  being 
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implemented  into  the  algorithm  and  is  briefly  described  below.  The  discussion  follows  closely  the 
discussion  presented  in  Hoteit  et  al.  [19]  and  is  presented  here  for  completeness.  For  further  details,  the 
reader  should  consult  the  references. 

Let  Ka  denote  an  arbitrary  element  and  K„  i  =1,2,3  the  neighbors  of  the  element.  Further,  let  i  =  0, 
1,2,3  denote  position  vectors  of  the  centroids  of  the  elements,  and  m,  denote  the  position  vectors  of  the 
midpoints  of  the  edges  of  element  Ka.  This  is  illustrated  in  Fig.  2. 


Fig.  2  —  Triangle  nomenclature  for  slope  limiting  procedure  of  Cockburn  and  Shu 


For  any  edge  ml5  we  can  write 


™\-b0  =al(bl-b0)  +  a2(b2-b0),  (29) 

where  (ai,  Ob)  are  constants  that  depend  only  on  the  geometry.  Additionally,  for  any  linear  function  «/„  we 
can  write 

uh  (mi  )~uh(b0)  =  04  (uh  (6, )  -  uh  (h0 ))  +  a2  (uh  (b2 )  -  uh  ( b0 )) .  (30) 

The  element  average  uh  is  the  value  of  the  function  «/,  evaluated  at  the  element  centroid.  Additionally, 
we  can  define 


uh  ( mhK0  )  =  uh  (m, )  -  u (K0 ), 


(31) 


which,  by  Eq.  (30),  can  be  written  as 


uh (m, , Kq )  =  ct! (uK]  -uKq)  +  a2 {uKi_  -  uKq  )  -  Am (m,  ,K0). 


(32) 


With  the  above  definitions,  the  slope  limiting  proceeds  as  described  below.  First  compute  the  quantity 
A,  for  each  edge  of  the  element  Ka  as 
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A i  =  M (uh (mi,K0 ),  v(A u (ml , K0 ))), 


(33) 


where  Mis  a  minmod  function  defined  by 


M(a]  ,a2) 


f 5  min  |  ctj  |  if  s  =  sign(a\ )  =  sign(a2 ) 
1 0  otherwise 


(34) 


and  v  is  a  constant  greater  than  one.  Cockbum  and  Shu  use  v  =  1.5.  With  A,  computed  by  Eq.  (33),  there 
are  two  possibilities: 

3 

1 .  If  7A:=0,  then 

i=i 


3 


An uh  =uKq  +  Z  a ih(x,y), 
i= 1 


(35) 


where  AYluh  is  the  slope-limited  value  of  Uh  and  §(x,y)  are  the  basis  functions. 

3 

2- If  I  A(  72  0 ,  compute  the  following: 


i=i 


3  3 

pos  =  Z  max(0,Af)  neg=  Z  max(0,-Af) 

i=l  i=l 

f  \  I  \ 


0  =  min 


neg 


0  =  min 


v  pos; 


pos 


V  neg; 


A,-  =  0+max(O,A;)  -0  max(0,-A;). 


Then  set 


3  „ 

KUuh  =  uKq  +  Z  A i§i(x,y). 


(36) 


(37) 


The  limiting  operator  locally  conserves  mass  and  guarantees  that  the  reconstructed  gradient  of  Anw,,  is 
not  greater  than  that  of  ip. 

5.  IMPLEMENTATION  INTO  QUODDY5 

The  details  of  the  numerical  implementation  of  the  DG  algorithm  are  described  below. 


12 


Guillot  et  al. 


5.1  Numerical  Algorithm 

QUODDY  contains  a  subroutine  called  ELEVATIONQ5  that  computes  the  depth  /fusing  the  GWCE 
and  a  continuous  Galerkin  finite-element  approximation.  The  numerical  implementation  of  the  DG 
method  replaces  ELEVATIONQ5  with  a  subroutine  ELEVATIONQ5DG,  which  serves  as  the  driver 
program  to  compute  the  depth  using  the  DG  method.  The  subroutine  calls  additional  subroutines  to 
complete  the  computations  described  in  Section  4. 

5.2  Projection  to  and  from  Nodes 

The  input  for  QUODDY5  provides  values  of  the  dependent  variables  at  the  nodes.  Elowever,  the  DG 
method  1  requires  the  degrees  of  freedom  to  be  defined  in  terms  of  the  element  average  and  the  derivative 
with  respect  to  x  and  y.  Additionally,  once  the  depth  is  computed  in  terms  of  the  DG  degrees  of  freedom, 
they  must  be  projected  back  to  the  nodes  so  that  they  can  be  used  to  compute  the  velocities  using  the 
continuous  Galerkin  method. 

The  projection  from  the  nodes  is  found  by  solving  the  3 

1  (x]  -xj)  O,  -yf) 

1  (x2~xj)  (y2  -vj) 

1  (x3  -xj)  (y3-yj) 


x  3  matrix. 


dy 


HJ0 

dHi 

Hi 

dx 

= 

h2 

(38) 

dHj 

[h3\ 

for  each  element.  Similar  projections  are  used  to  compute  the  DG  degrees  of  freedom  for  U  and  V. 

After  computing  the  depth  using  the  DG  method,  it  is  necessary  to  project  the  DG  degrees  of  freedom 
back  onto  the  nodes.  This  is  not  as  straightforward  as  projecting  from  the  nodes,  since  each  node  can  be 
connected  to  any  number  of  elements.  The  projection  back  onto  the  nodes  is  found  by  computing  the 
value  of  the  depth  at  the  node  using  the  DG  degrees  of  freedom  for  each  element  attached  to  the  node,  and 
then  dividing  by  the  total  number  of  elements  attached  to  the  node.  The  algorithm  proceeds  as  follows. 
Let  ek  denote  an  element  attached  to  node  n,  and  H„(ek)  the  elevation  computed  at  node  n  using  the  basis 
functions  on  ek.  If  a  total  of  K  elements  are  attached  to  node  n,  then  the  average  value  H„  at  the  node  is 
given  by 


Hn=^THn(ek). 

K  k=\ 


(39) 


5.3  Input  Modification 

The  DG  implementation  into  QUODDY  changes  the  input  minimally.  The  idea  was  to  make  the 
implementation  as  transparent  to  the  user  as  possible.  Only  one  file,  the  .inq  file,  has  been  changed.  A 
typical  .inq  is  shown  in  Fig.  3.  Only  one  line  has  been  added  to  the  file.  The  last  line  is  a  flag  that  allows 
the  user  to  choose  “GWCE”  to  use  the  original  formulation  or  “DGSL”  to  choose  the  DG  method.  Within 
the  code  itself,  anyplace  where  modification  to  existing  code  has  been  made,  the  lines  mjg_begin  and 
mjg  end  have  been  placed  between  the  modified  code.  Finally,  a  file  called  quoddy_5.1.1_DG_SUBS.f 
that  contains  the  DG  subroutines  has  been  added  to  the  existing  suite  of  subroutines.  This  file  must  be 
added  to  the  makefile  and  compiled  with  the  existing  code. 
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( Comment : ) 

bankl50  barotropic  tidal  solution 
(Case  name: ) 

. ./banklSO/banklSO 
(Coordinates:  CARTESIAN/SPHERICAL) 

CARTESIAN 

(Boundary  element  incidence  list:) 

. . /bankl50/bankl50_elev .bel 
(Initial  condition  file:) 

COLD-START  01  01  1996  0.0 
(Echo  file : ) 
bankl50_elev.CT.echo 
(Simulation  parameters:) 

SI  UNITS  (units) 

1.00  1.00  1.00  (x,  y,  and  z  scaling  factor) 

43.500  (degree  latitude] 

10.0  (minimum  depth) 

01  01  1996  536544.0  [end  date  (d  m  y)  and  time  (sec)  of  simulation) 
15.65625  [time  step  (seconds)] 

21  [number  of  vertical  nodes) 

GWCE  IGWCE  or  DGSL] 


Fig.  3  — Typical  .inq  file 


6.  VALIDATION 

The  implementation  is  validated  in  the  test  cases  described  below. 


6.1  Maximum  Courant  Number  Verification 

Before  choosing  a  test  case,  the  courant  number  restriction  was  verified  for  a  trivial  case  in  which  the 
surface  elevation  is  initially  zero  and  zero  elevation  boundary  conditions  are  imposed.  The  elevation 
throughout  the  domain  should  remain  zero  for  all  t  >  0.  The  Gulf  of  Maine  mesh  is  used  in  this  case.  The 
maximum  courant  number  for  each  element  is  shown  in  Fig.  4  along  with  the  instability  in  the  solution 
after  40  time  steps,  using  a  time  step  At  =  25  s.  If  the  time  step  is  reduced  so  that  the  maximum  courant 
number  in  the  mesh  is  below  unity,  the  solution  remains  stable. 

The  implementation  was  tested  on  a  case  provided  by  the  developers  of  QUODDY.  The  test  case  is  the 
Gulf  of  Maine  with  periodic  tidal  (elevation)  boundary  conditions  applied  on  all  boundaries.  The  mesh 
and  bathymetry  are  shown  in  Fig.  5.  The  depth  is  in  meters. 

6.2  Gulf  of  Maine  with  Elevation  Boundary  Conditions 

The  elevation  was  computed  using  both  the  DG  and  GWCE  methods.  The  tidal  boundary  conditions 
were  ramped  up  over  three  tidal  periods  and  a  total  of  1 0  tidal  periods  were  computed.  The  slopes  of  the 
total  depth  are  held  constant  at  their  initial  values  in  the  computations.  A  less  restrictive  slope  limiter 
developed  by  Cockbum  and  Shu  is  currently  being  implemented  into  the  algorithm. 

The  elevations  computed  using  the  DG  and  GWCE  methods  are  shown  at  3  days  along  with  the 
differences  in  elevation  in  Fig.  6.  The  computed  elevation  solutions  are  very  similar  and  generally  within 
1%  to  2%.  The  corresponding  velocity  vectors  are  shown  in  Fig.  7.  The  velocity  vectors  are  also 
quantitatively  very  similar. 
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6.3  Discussion  of  Different  Methods  Investigated 

Two  methods  of  computing  the  degrees  of  freedom  per  element  were  investigated.  In  the  first,  the 
element  degrees  of  freedom  were  the  element  average  and  gradient  and  in  the  second,  they  were  the 
values  at  the  element  vertices.  It  was  found  that  each  method  produced  similar  results,  although,  in  the 
second  method,  computing  elevation  boundary  conditions  was  facilitated  by  simply  setting  the  value  at 
the  given  boundary  node  to  the  value  given  by  the  prescribed  boundary  condition.  Additionally,  the 
element  vertices  method  produced  a  simpler  mass  matrix  that  reduced  the  computational  effort  required. 
In  addition  to  the  methods  investigated,  Cockbum  and  Shu  [15]  used  the  element  edge  midpoints  as  the 
degrees  of  freedom.  This  has  the  advantage  of  producing  a  diagonal  mass  matrix  for  each  element. 
However,  in  any  case,  the  mass  matrix  is  a  3  x  3  matrix  that  is  easily  inverted  by  hand,  so  this  is  not  a 
great  advantage,  and  Hoteit  et  al.  [19]  noted  that  using  the  midpoints  of  the  edges  produces  “excessive 
smearing”  of  the  solution. 


Courant  Number 


Fig.  4  —  Courant  number  (left)  and  instability  in  elevation  solution  (right)  after  40  time  steps,  At  =  25 


Fig.  5  — Gulf  of  Maine  mesh  and  bathymetry  (left)  plan  view  and  (right)  3-D  view 
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Difference  in  Solutions 


Fig.  6  —  Elevation  solution  in  meters  at  t  =  3  days  using  DG 
and  GWCE  methods,  and  difference  in  solutions 
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Difference  in  V-Velocity 


Fig.  7  —  Velocity  solution  in  m/s  at  t  =  3.0  days  using  DG  and 
GWCE  methods,  and  difference  in  solutions 


Two  different  second-order  Runge-Kutta  methods  were  also  investigated.  It  was  found  that  Henn’s 
method,  described  by  Eq.  (17),  produced  somewhat  more  stable  results  than  the  method  described  by  Eq. 
(16). 
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While  both  methods  produced  unstable  results  if  the  slopes  were  not  restricted,  the  instability  using 
Henn’s  method  grew  at  a  slower  rate  than  it  did  using  the  other  method. 

7.  CONCLUSIONS 

A  discontinuous  Galerkin  method  was  implemented  into  QUODDY5  to  compute  the  depth  using  the 
primitive  conservation  of  mass  equation,  replacing  the  GWCE  formulation  in  the  code.  Several  variations 
were  investigated  and  the  advantages  and  disadvantages  of  each  were  discussed. 

The  explicit  formulation  of  the  DG  method  limits  the  time  step  because  of  a  CFL  condition  as 
compared  to  the  semi-implicit  formulation  of  the  GWCE  method.  This  time-step  restriction  is  problem 
and  mesh  dependent,  but  for  the  Gulf  of  Maine  case,  the  maximum  allowable  time  step  of  the  DG  method 
was  approximately  1/10  of  the  time  step  using  the  GWCE  method,  resulting  in  longer  run  times  for  the 
same  problem. 

The  DG  method  requires  a  slope  limiter  to  ensure  stability  of  the  solution.  The  TVB  slope  limiter  of 
Cockbum  and  Shu  is  currently  being  implemented  into  the  code. 

8.  RECOMMENDATIONS 

Future  work  in  this  area  should  include  developing  a  fully  DG  version  of  QUODDY.  This  will 
facilitate  a  comparision  of  the  DG  method  using  the  primitive  conservation  of  mass  equation  with  the 
continuous  Galerkin  method  using  the  GWCE  method.  Since  strengths  of  the  DG  method  include  ease  of 
parallelization  and  h-p  adaptivity,  longer  term  goals  should  include  developing  a  parallel  version  of  the 
method  including  dynamic  mesh  adaptivity. 
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